R on statistikale ja andmete visualiseerimisele suunatud vabavaraline programmeerimiskeel. R loodi 90ndate alguses S-i baasile ning on tasapisi muutunud üheks enimkasutatavaks statistika tarkvaraks.
R Passes SPSS in Scholarly Use, Stata Growing Rapidly (2014)
R-il on olulisi eeliseid mitmete teiste statistikatarkvarade ees:
Sellega seoses on sel ka miinuseid:
Sissejuhatus:
Jaotus ja keskmised:
Usaldusvahemikud:
Keskmiste võrdlus:
Visualiseerimine:
Lineaarne regressioon:
Logistiline regressioon:
Huvitavat:
R-i kasutamiseks piisab vaid R-i “mootori” installeerimisest, kuid mugavam on kasutada graafilist kasutajaliidest, näiteks R-Studiot.
Kuna R on vabavaraline tarkvara, siis on mitmete R-i funktsioonide kasutamiseks vajalik need esmalt alla laadida. See toimub R-i keskkonnas install.packages(...) abil (nt install.packages("tidyverse")), mispuhul pakett on funktsioonide kogum. Vajaminevate pakettide peale hetkel mõtlema ei pea, need tulevad jooksvalt praktikumide käigus.
R-i pakettidel on enamasti küllatki hästi lahti kirjutatud help sektsioon, mis on leitav help tabi alt (enamasti paremal all nurgas). Help tabi otsingusse saab kirjutada otsitud funktsiooni nime, misjärel ta kuvab ta funktsiooni kirjelduse ning argumentide väärtused. Enamasti on lisatud mõned näited funktsiooni kasutamisvõimalustest. Võimalik on veel:
?read.csv()) või help(...), kuhu läheb argumendiks funktsiooni nimi;Lisaks R-is olemasolevale kirjeldusele on mitmed paketid avaldatud teadusajakirjades, eeskätt ajakirjas Journal of Statistical Software. Samuti on enamus tüüp-probleemid lahendatud Stack Overflow lehel. Mõlemal juhul peaks google viima sobiva lahenduseni.
R-i põhiline töövahend on objekt, mis võib sisaldada nii funktsioone kui andmeid. Objekte saab ise luua <- või = märkide abil (<- on mõnevõrra kindlam meetod ja mõned soovitavad kasutada seda. Üldjuhul, sh selle aine raames, piisab ka = kasutamisest).
x = 5
x
[1] 5
Kui soovime siduda rohkem andmeid ühe objektiga, on võimalik luua vektor (jada / ühedimensionaalne maatriks).
x <- c(5, 12.3, 3:10)
x
[1] 5.0 12.3 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
Objekt võib saada sisendi ka teisest objektidest. Samuti võib teha objektidega tehteid.
y = x - 5
y
[1] 0.0 7.3 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0
Objektiks võib olla ka funktsioon.
lahuta.viis = function(x){
tmp = x - 5
return(tmp)
}
y = lahuta.viis(x)
y
[1] 0.0 7.3 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0
Numbrid:
x = 1
str(x)
num 1
Tähemärgid (numbrid muudetakse tähemärkideks):
x = c("a", "b", 1, 2)
str(x)
chr [1:4] "a" "b" "1" "2"
Faktorid (identsed elemendid saavad sama koodi):
x = as.factor(c("a", "b", "1", "2", "a", "2"))
str(x)
Factor w/ 4 levels "1","2","a","b": 3 4 1 2 3 2
Kuupäevad:
x = as.Date("2017-09-01")
str(x)
Date[1:1], format: "2017-09-01"
Kuupäevade puhul on oluline jälgida formaati. Vaikimisi on selleks yyyy-mm-dd, kuid kui sisse loetavad andmed on teises formaadis, tuleks see as.Date(...) käskluses täpsustada format või origin argumentidega. Formaatidest pikemalt: Date Values.
as.Date("01-09-2017", format="%d-%m-%Y")
[1] "2017-09-01"
as.Date(17410, origin="1970-01-01")
[1] "2017-09-01"
Seni vaadatud jadad on olnud kõik ühedimensionaalsed, ent andmed võivad olla ka mitmedimensionaalsed. Enimlevinud mitmedimensionaalsed formaadid on andmetabelid ja maatriksid, millest peamiselt keskendume esimesele. Maatriks on kahedimensionaalne vektor, mis on justkui mitu vektorit kokku liidetuna, kuid kõik vektorid peavad olema sama liiki andmetega.
matrix(1:9, nrow=3, byrow = T)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
matrix(letters[1:9], nrow=3, byrow = T)
[,1] [,2] [,3]
[1,] "a" "b" "c"
[2,] "d" "e" "f"
[3,] "g" "h" "i"
Mitu samapikkust vektorit saab liita kokku üheks ka andmetabeliks ehk data frameiks (data.frame(...)). Data frame’i erinevus seisneb selles, et veerud võivad sisaldada eritüüpi andmeid (nt üks veerg on numbriline, teine sisaldab faktoreid jne).
x = c(1:4)
y = c("Audi", "DAF", "Ford", "GAZ")
z = x-5
df = data.frame(tunnus.1 = x,
tunnus.2 = y,
tunnus.3 = z)
df
tunnus.1 tunnus.2 tunnus.3
1 1 Audi -4
2 2 DAF -3
3 3 Ford -2
4 4 GAZ -1
Andmetabelite tunnuseid (ehk veerge) saab kätte ka vektorina kasutades $ märki:
df$tunnus.2
[1] Audi DAF Ford GAZ
Levels: Audi DAF Ford GAZ
Võimalik on kasutada ka indeksi (järjekorra) numbrit kandilistes sulgudes. Andmetabelite puhul on esimene number rea number ning koma järel veeru number. Vektori puhul piisab vaid rea numbrist.
df[,2]
[1] Audi DAF Ford GAZ
Levels: Audi DAF Ford GAZ
df["tunnus.2"]
tunnus.2
1 Audi
2 DAF
3 Ford
4 GAZ
df[1,]["tunnus.2"]
tunnus.2
1 Audi
df[1,2]
[1] Audi
Levels: Audi DAF Ford GAZ
df$tunnus.2[2]
[1] DAF
Levels: Audi DAF Ford GAZ
Kõiki erinevaid objekte saab koondada listidesse (list(...)). List on hea koondamaks erinevate pikkustega või erineva iseloomuga, kuid samateemalisi andmeid. Tihti võib liste kohata erinevate funktsioonide väljundites.
ls = list(df)
ls[[2]] = c(1:10)
ls[[3]] = y[2:4]
ls
[[1]]
tunnus.1 tunnus.2 tunnus.3
1 1 Audi -4
2 2 DAF -3
3 3 Ford -2
4 4 GAZ -1
[[2]]
[1] 1 2 3 4 5 6 7 8 9 10
[[3]]
[1] "DAF" "Ford" "GAZ"
Kõikide objektide puhul:
str(...) - objekti struktuurtypeof(...) - storage modeis.numeric(...) - kas on…is.character(...)is.factor(...)is.function(...)is.logical(...)…
Vektori puhul:
length(...) - vektori pikkusData frame’i puhul:
dim(...) - data frame’i dimensioonid (ridu, veerge)nrow(...) - ridade arvncol(...) - veergude arvFormaadi muutmiseks:
as.vector(...) - vektoras.matrix(...) - maatriksas.Date(...) - kuupäevas.data.frame(...) - data frameas.list(...) - listR suudab lugeda sisse enamikke levinud andmefailide formaate, kuid kõige sagedasemalt kasutatakse .csv (comma separated values). Andmestike sisse lugemiseks on erinevaid funktsioone, kuid .csv-de puhul piisab read.csv(...). read.csv(...) võimaldab muu hulgas määratleda, kuidas on eraldatud erinevad veerud (sep; koma, semikooloni vm abil) ja mida kasutatakse komakohtade eraldamiseks (dec; punkt, koma vm). See võib olla oluline, kuna erinevad programmid erinevatel OS-idel võivad kasutada erinevaid standardeid.
url = "https://martenveskimae.github.io/stat-mod/party_dataset.csv"
df = read.csv(url)
Andmestiku vaatamiseks võib küsida tervet andmestikku, aga võib ka küsida esimesi või viimaseid ridu. Selleks saab kasutada nii indeksi kui ka eraldi funktsioone head(...) ja tail(...).
df[1:6, 1:4]
party support support_emor last
1 ref 0.36 NA 0.278
2 ref 0.34 NA 0.278
3 ref 0.45 NA 0.278
4 ref 0.43 NA 0.278
5 ref 0.43 NA 0.278
6 ref 0.44 NA 0.278
head(df[,1:4])
party support support_emor last
1 ref 0.36 NA 0.278
2 ref 0.34 NA 0.278
3 ref 0.45 NA 0.278
4 ref 0.43 NA 0.278
5 ref 0.43 NA 0.278
6 ref 0.44 NA 0.278
tail(df[,1:4])
party support support_emor last
551 sde NA 0.15 0.152
552 irl NA 0.06 0.137
553 kesk NA 0.23 0.248
554 ref NA 0.26 0.277
555 vaba NA 0.07 0.087
556 ekre NA 0.16 0.081
Et näha, millised tunnused andmestikus on, võib küsida tunnuste pealkirju colnames(...) funktsiooiga.
colnames(df)
[1] "party" "support"
[3] "support_emor" "last"
[5] "year" "month"
[7] "gov" "month_cycle"
[9] "cycle_constant" "cycle"
[11] "inflats_m_basevalue" "inflats_m"
[13] "unemp" "gdpgr"
[15] "gdpgr_yearly" "skp_yearly"
[17] "support_lag1" "bruto"
[19] "bruto_yearly" "inflats_prevyear_base"
[21] "alko_suits_base" "alko_suits_prevyear"
[23] "eluase_base" "eluase_prevyear"
[25] "irl_meedias" "kesk_meedias"
[27] "ekre_meedias" "ref_meedias"
[29] "sde_meedias" "rohel_meedias"
[31] "date"
Kohalikult kettalt andmestike laadimiseks tuleb märkida ära andmestiku asukoht kettal. Selle lihtsustamiseks on otstarbekas märkida ära töökeskonna asukoht.
Olemasoleva asukoha nägemiseks on funktsioon getwd().
getwd()
[1] "/Users/martenveskimae/GitHub/stat-mod"
Kui tahta seda muuta, saab kasutada setwd(...) funktsiooni, kus argumendiks on soovitud uus asukoht.
Kokkuvõtvaid andmed nii andmetabelite kui teiste andmeformaatide kohta näeb funktsiooniga summary(...).
summary(df)
party support support_emor last
ekre: 30 Min. :0.040 Min. :0.0500 Min. :0.0810
irl :124 1st Qu.:0.130 1st Qu.:0.1000 1st Qu.:0.1520
kesk:124 Median :0.210 Median :0.1600 Median :0.2100
ref :124 Mean :0.204 Mean :0.1601 Mean :0.2014
sde :124 3rd Qu.:0.260 3rd Qu.:0.2200 3rd Qu.:0.2610
vaba: 30 Max. :0.450 Max. :0.2800 Max. :0.2900
NA's :34 NA's :418 NA's :12
year month gov month_cycle
Min. :2007 Min. : 1.000 Min. :0.0000 Min. : 1
1st Qu.:2010 1st Qu.: 3.000 1st Qu.:0.0000 1st Qu.:11
Median :2012 Median : 6.000 Median :1.0000 Median :20
Mean :2012 Mean : 6.378 Mean :0.5665 Mean :22
3rd Qu.:2015 3rd Qu.: 9.000 3rd Qu.:1.0000 3rd Qu.:33
Max. :2017 Max. :12.000 Max. :1.0000 Max. :49
NA's :12
cycle_constant cycle inflats_m_basevalue inflats_m
Min. :49 Min. :0.02041 Min. :108.5 Min. :-0.90701
1st Qu.:49 1st Qu.:0.22449 1st Qu.:124.7 1st Qu.:-0.03494
Median :49 Median :0.40816 Median :138.0 Median : 0.28775
Mean :49 Mean :0.44898 Mean :134.2 Mean : 0.30774
3rd Qu.:49 3rd Qu.:0.67347 3rd Qu.:143.9 3rd Qu.: 0.64917
Max. :49 Max. :1.00000 Max. :146.3 Max. : 2.08868
NA's :12 NA's :12 NA's :72 NA's :180
unemp gdpgr gdpgr_yearly
Min. : 4.000 Min. :-10.80000 Min. :-18.4000
1st Qu.: 6.500 1st Qu.: -0.30000 1st Qu.: -1.7000
Median : 8.000 Median : 0.70000 Median : 2.3000
Mean : 9.141 Mean : -0.01209 Mean : 0.7391
3rd Qu.:11.300 3rd Qu.: 1.30000 3rd Qu.: 5.6000
Max. :19.500 Max. : 3.70000 Max. : 11.6000
NA's :78 NA's :192 NA's :162
skp_yearly support_lag1 bruto bruto_yearly
Min. :-19.3000 Min. :0.0700 Min. : 659.7 Min. :-6.500
1st Qu.: 0.3000 1st Qu.:0.1500 1st Qu.: 792.3 1st Qu.: 3.900
Median : 1.9000 Median :0.2300 Median : 857.1 Median : 6.300
Mean : 0.8247 Mean :0.2231 Mean : 898.4 Mean : 6.707
3rd Qu.: 5.2000 3rd Qu.:0.2800 3rd Qu.:1010.0 3rd Qu.: 8.500
Max. : 10.5000 Max. :0.4500 Max. :1163.0 Max. :21.166
NA's :78 NA's :184 NA's :84 NA's :192
inflats_prevyear_base alko_suits_base alko_suits_prevyear eluase_base
Min. :102.7 Min. :106.8 Min. :102.1 Min. :121.7
1st Qu.:123.2 1st Qu.:138.1 1st Qu.:115.3 1st Qu.:149.6
Median :133.1 Median :149.3 Median :141.2 Median :168.3
Mean :130.6 Mean :148.9 Mean :139.5 Mean :168.3
3rd Qu.:143.3 3rd Qu.:167.3 3rd Qu.:158.3 3rd Qu.:191.4
Max. :145.7 Max. :182.4 Max. :174.9 Max. :199.0
NA's :60 NA's :174 NA's :174 NA's :126
eluase_prevyear irl_meedias kesk_meedias ekre_meedias
Min. :106.1 Min. : 227.0 Min. : 312.0 Min. : 8.0
1st Qu.:144.7 1st Qu.: 438.5 1st Qu.: 559.8 1st Qu.: 29.0
Median :158.0 Median : 559.0 Median : 626.5 Median : 51.5
Mean :159.5 Mean : 610.4 Mean : 717.0 Mean :115.7
3rd Qu.:183.2 3rd Qu.: 730.2 3rd Qu.: 811.5 3rd Qu.:111.0
Max. :199.0 Max. :1597.0 Max. :2414.0 Max. :673.0
NA's :126 NA's :316 NA's :316 NA's :316
ref_meedias sde_meedias rohel_meedias date
Min. : 343.0 Min. : 245.0 Min. : 8.00 2015-01-01: 6
1st Qu.: 537.8 1st Qu.: 367.0 1st Qu.: 26.75 2015-02-01: 6
Median : 658.5 Median : 444.5 Median : 49.00 2015-03-01: 6
Mean : 747.7 Mean : 510.6 Mean : 85.40 2015-04-01: 6
3rd Qu.: 929.0 3rd Qu.: 585.8 3rd Qu.:103.00 2015-05-01: 6
Max. :1836.0 Max. :1480.0 Max. :462.00 2015-06-01: 6
NA's :316 NA's :316 NA's :316 (Other) :520
Faktortunnuste puhul võib pakkuda huvi, millised erinevad kategooriad on andmestikus esindatud. Näiteks võime vaadata, millised erakonnad andmestikus leiduvad.
levels(df$party)
[1] "ekre" "irl" "kesk" "ref" "sde" "vaba"
Kui soovida täpsemat ülevaadet, näiteks, mitu andmepunkti on iga erakonnaga seotud, siis selleks on hea kasutada paketti tidyverse. Tidyverse hõlmab endas mitmeid erinevaid, kuid omavahel ühilduvaid pakette. Üks neist, dplyr, pakub mitmeid funktsioone andmete mugavaks töötlemiseks (eeskätt filtreerimine, grupeerimine, tunnuste muutmine jne). Teine oluline eelis tuleneb paketist magrittr, mis võimaldab käsklusi sisestada jadana (ing.k pipe). See kiirendab tööprotsessi ja muudab R-i koodi kergemini jälgitavaks. Pipe’ide kasutamiseks tuleb kahe funktsiooni vahel lisada %>% sümbolid. Enne näite juurde liikumist on vajalik tidyverse aktiveerida, milleks on funktsioon library(...).
library(tidyverse)
Pipe’i näide võib olla järgnev:
c(1:10) %>% lahuta.viis()
[1] -4 -3 -2 -1 0 1 2 3 4 5
Järgnevalt tutvume, mitu rida on seotud üksikute erakondadega (sagedustabel).
df %>%
group_by(party) %>%
summarize(n = n())
# A tibble: 6 x 2
party n
<fctr> <int>
1 ekre 30
2 irl 124
3 kesk 124
4 ref 124
5 sde 124
6 vaba 30
dplyr kohta leiab rohkem informatsiooni siit: Grammar of Data Manipulation • dplyr
Kõige kiirem viis andmete visualiseerimiseks on kasutada funktsioone plot(...) või hist(...), kus tuleb määratleda telje/telgede väärtused.
plot(df$year, df$support)
hist(df$support)
Paraku ei ole antud funktsioonid väga mugavad tidyverse’is kasutamiseks ning kujunduse muutmine on piiratud võimalustega. Kõige enim levinud viis andmete visualiseerimiseks R-i keskkonnas on paketiga ggplot2, mis on samuti tidyverse’i osa.
df %>%
filter(party=="ref") %>%
ggplot() +
geom_point(aes(year, support))
df %>%
filter(party=="ref") %>%
ggplot() +
geom_histogram(aes(support))
Võrdlusoperaatorid
Filtreerimisel või muude võrdluste tegemisel kasutatakse teatud tingmärke, mida nimetatakse võrldusoperaatoriteks. R-is on need järnevad:
| Operaator | Tulemus |
|---|---|
== |
Võrdne |
!= |
Mittevõrdne |
> |
Suurem |
< |
Väiksem |
>= |
Suurem-võrdne |
<= |
Väiksem-võrdne |
%in% |
Sisaldab |
Loogilised operaatorid
Olukordades, kus soovime mitut võrdlust/tingimust kombineerida, tuleb erinevad komponendid siduda omavahel loogiliste operaatoritega:
| Operaator | Tulemus |
|---|---|
& |
Loogiline JA üksikute vektori komponentide kohta (väljund on TRUE/FALSE kõigi vektori komponentide kohta) |
&& |
Loogiline JA kõigi vektori komponentide kohta (väljund on ainult üks TRUE/FALSE kogu vektori kohta) |
| |
Loogiline VÕI üksikute vektori komponentide kohta (väljund on TRUE/FALSE kõigi vektori komponentide kohta) |
|| |
Loogiline VÕI kõigi vektori komponentide kohta (väljund on ainult üks TRUE/FALSE kogu vektori kohta) |
! |
Loogiline mittesamaväärsus |
Näiteid: R - Operators
Kujunduse kiireks muutmiseks on erinevad valmislahendused, algusega theme_.
df %>%
filter(party=="ref") %>%
ggplot() +
geom_point(aes(year, support)) +
theme_bw()
ggplot2 kohta leiab rohkem informatsiooni siit: Elegant Data Visualisations Using the Grammar of Graphics
library(tidyverse)
library(e1071)
library(reshape2)
df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
Vaatame jaotuste järsakust ja sümmeetrilisust visuaalselt histogrammi abil. geom_histogram(...) annab vaikimisi sagedusjaotuse, kuid on võimalik valida ka protsente või tõenäosustihedust.
df %>%
ggplot() +
geom_histogram(aes(support, ..count../sum(..count..))) +
geom_line(aes(support, ..density../100), stat = "density") +
scale_y_continuous(labels=scales::percent) +
scale_x_continuous(labels=scales::percent) +
labs(y="protsent") +
theme_bw()
ggplot võimaldab luua eraldi joonised erinevate gruppide jaoks kasutades facet_wrap(...) käsklust.
df %>%
ggplot() +
geom_histogram(aes(support, ..count../sum(..count..))) +
geom_line(aes(support, ..density../100), stat = "density") +
facet_wrap(~party) +
scale_y_continuous(labels=scales::percent) +
scale_x_continuous(labels=scales::percent) +
labs(y="protsent") +
theme_bw()
Silmas tuleb pidada, et sel puhul kalkuleerib ggplot tulpade laiused ja hulga üle kõikide kuvatud jooniste, mistõttu erinevad need üksikult kuvatud joonistest. Tulpade hulga ja laiuste kontrollimiseks on argumendid bins ja binwidth, millest esimene määrab tulpade hulga ja teine laiused. Neid käsklusi ei saa korraga kasutada (binwidth kirjutab binsi üle).
df %>%
ggplot() +
geom_histogram(aes(support, ..count../sum(..count..)), bins=10, binwidth=0.03) +
geom_line(aes(support, ..density../100), stat = "density") +
facet_wrap(~party) +
scale_y_continuous(labels=scales::percent) +
scale_x_continuous(labels=scales::percent) +
labs(y="protsent") +
theme_bw()
Vaatame järsakust ja sümmeetrilisust ka numbriliselt nn järsakuse- ja assümmeetriakordaja abil. Seda võimaldab pakett e1071, kust leiame käsklused skewness(...) ja kurtosis(...) (rakendatavad valemid on leitavad helpist). Kuna soovime rakendada antud käsklusi kõigi gruppide (erakondade) jaoks eraldi, siis tuleks andmestik esmalt filtreerida grupeerimistunnuse alusel (party) ning alles seejärel rakendada antud käsklusi.
df %>%
group_by(party) %>%
summarise(skewness = skewness(support, na.rm=T),
kurtosis = kurtosis(support, na.rm=T))
# A tibble: 6 x 3
party skewness kurtosis
<fctr> <dbl> <dbl>
1 ekre -0.23712096 -0.8569726
2 irl -0.61654227 -0.2531924
3 kesk 0.29345645 -0.2409381
4 ref 0.30822221 -0.5851121
5 sde 0.44024113 -1.2410302
6 vaba 0.01834276 -1.3350304
Alternatiiv dplyr meetodile oleks for loop. Loopimine on küll võrdlemisi lihtne, kuid vahel mitte kõige efektiivsem viis:
parties = as.character(unique(df$party))
s = c()
k = c()
for(i in 1:length(parties)){
s[length(s)+1] = skewness(df$support[df$party==parties[i]], na.rm=T)
k[length(k)+1] = kurtosis(df$support[df$party==parties[i]], na.rm=T)
}
names(s) = parties
names(k) = parties
s;k
ref kesk irl sde vaba ekre
0.30822221 0.29345645 -0.61654227 0.44024113 0.01834276 -0.23712096
ref kesk irl sde vaba ekre
-0.5851121 -0.2409381 -0.2531924 -1.2410302 -1.3350304 -0.8569726
Sama tulemuse saab efektiivsemalt kasutades apply(...) käsklust (antud juhul sapply(...). Rohkem applyst: Tutorial on the R Apply Family), mis jäljendab eelneva meetodi loogikat.
s = sapply(as.character(unique(df$party)), function(x) skewness(df$support[df$party==x], na.rm=T))
k = sapply(as.character(unique(df$party)), function(x) kurtosis(df$support[df$party==x], na.rm=T))
s;k
ref kesk irl sde vaba ekre
0.30822221 0.29345645 -0.61654227 0.44024113 0.01834276 -0.23712096
ref kesk irl sde vaba ekre
-0.5851121 -0.2409381 -0.2531924 -1.2410302 -1.3350304 -0.8569726
Moodi, ehk kõige sagedamini esineva väärtuse leidmiseks eraldi käsklus puudub, ent seda on lihtne leida sagedustabelist.
df %>%
group_by(party, support) %>%
summarise(n = n()) %>%
na.omit() %>%
group_by(party) %>%
summarise(mode = support[n==max(n)][1])
# A tibble: 6 x 2
party mode
<fctr> <dbl>
1 ekre 0.10
2 irl 0.16
3 kesk 0.26
4 ref 0.25
5 sde 0.12
6 vaba 0.11
dplyri saab kasutada ka keskmise ja mediaani leidmiseks.
df %>%
group_by(party) %>%
summarise(mean = mean(support,na.rm=T),
median = median(support,na.rm=T))
# A tibble: 6 x 3
party mean median
<fctr> <dbl> <dbl>
1 ekre 0.10952381 0.10
2 irl 0.13875000 0.15
3 kesk 0.25858333 0.26
4 ref 0.28716667 0.29
5 sde 0.16633333 0.14
6 vaba 0.09904762 0.10
Viimaks võib kanda leitud väärtused histogrammile, et neid visuaalselt kõrvutada. Kuigi ggplot oskab lisada joonistele funktsioone, ei suuda ta praeguses versioonis kuvada eraldi funktsioone eraldi gruppidele (nt kasutades facet_wrap(~party)). Otstarbekas on sellisel juhul teha uus andmetabel keskmiste väärtustega. Andmetabeli mugavamaks kasutamiseks pöörame tabeli laiast formaadist pikka formaati, kasutades selleks reshape2 paketti ja melt(...) funktsiooni. Rohkem reshape2 ja laia/pika formaadi kohta: An Introduction to reshape2.
keskmised = df %>%
group_by(party, support) %>%
summarise(n = n()) %>%
na.omit() %>%
group_by(party) %>%
summarise(keskmine = mean(support),
mediaan = median(support),
mood = support[n==max(n)][1]) %>%
melt("party",c("keskmine", "mediaan", "mood"),"Keskmised")
df %>%
ggplot() +
geom_histogram(aes(support)) +
geom_line(aes(support, ..density..), stat = "density") +
geom_vline(data=keskmised,aes(xintercept=value, linetype=Keskmised, color=Keskmised)) +
facet_wrap(~party) +
scale_x_continuous(labels=scales::percent) +
theme_bw()
Lisaks mediaanile võivad huvipakkuvad olla ka teised protsentiilid (mediaan on 50 protsentiil), eeskätt 25 ja 75 (teisisõnu I ja III kvartiil). Veel võib aeg-ajalt kohta mõistet kvartiilhaare (interquartile range - IQR), mis tähistab vahemikku I ja III kvartiili vahel. R-is on protsentiilide leidmiseks käsklus quantile(...).
sapply(as.character(unique(df$party)), function(x) quantile(df$support[df$party==x],c(.025,.25,.5,.75,.975),na.rm=T))
ref kesk irl sde vaba ekre
2.5% 0.17975 0.2000 0.04975 0.08000 0.065 0.07
25% 0.23000 0.2375 0.12000 0.11750 0.080 0.10
50% 0.29000 0.2600 0.15000 0.14000 0.100 0.10
75% 0.33000 0.2800 0.16000 0.23000 0.120 0.13
97.5% 0.43000 0.3300 0.20025 0.28025 0.135 0.14
Eelnevalt toodud protsentiilid saab graafida boxploti abil, mis kuvab veel lisaks erindid (outliers).
df %>%
ggplot() +
geom_boxplot(aes(party,support)) +
scale_y_continuous(labels=scales::percent) +
theme_bw()
library(tidyverse)
library(reshape2)
df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
Eelmisest korrast teame, et keskmine toetusprotsent erakonniti erineb, aga erineb ka kõikumine keskmise ümber. Osade erakondade puhul on keskmine toetus stabiilselt ühe tüüpväärtuse ümber, aga teistel mitte. Selle hindamiseks kasutatakse tüüpiliselt variatsiooni ja standardhälvet, mis on R-is leitav käsklustega var(...) ja sd(...).
Variatsioon on eeldatav erinevus tõelise väärtuse ja keskmise vahel ruudus:
\(Var(X) = E[(X-\mu)^2]\)
ehk teisisõnu erinevuste keskmine:
\(Var(X) = \frac{ \sum(X-\mu)^2}{N-1}\)
X = 1:10
var(X)
[1] 9.166667
sum((X - mean(X))^2) / (length(X)-1)
[1] 9.166667
(Miks on jagaja \(N-1\) ja mitte \(N\)? Vahe tuleneb sellest, kas meil on andmed kogu populatsiooni kohta või valim populatsioonist, mis on alati väiksem kui N!)
Variatsiooni leidmiseks on erinevused võetud ruutu. Et jõuda tagasi algskaala juurde, tuleks variatsioonist võtta ruutjuur, mis ongi standardhälve:
\(\sigma=\sqrt{E[(X-\mu)^2]}\)
X = 1:10
sqrt(var(X))
[1] 3.02765
sd(X)
[1] 3.02765
Standardhälvet võib seega mõista kui keskmist erinevust keskmisest.
Vaatame erakondade keskmisi ja standardhälvet lähemalt:
keskmised = df %>%
group_by(party) %>%
summarise(keskmine = mean(support, na.rm=T),
s.ylemine = keskmine + sd(support, na.rm=T),
s.alumine = keskmine - sd(support, na.rm=T)) %>%
melt("party",c("keskmine", "s.ylemine", "s.alumine"), "Varieeruvus")
df %>%
ggplot() +
geom_line(aes(as.Date(date), support)) +
geom_hline(data=keskmised,aes(yintercept=value, linetype=Varieeruvus, color=Varieeruvus)) +
facet_wrap(~party) +
scale_y_continuous(labels=scales::percent) +
labs(x="") +
theme_bw()
Millest see kõikumine keskmine ümber tuleneb on meile teada. Erakondade toetus varieerub ajas erinevalt ja meie arvutatud keskmine näitab tüüpväärtust üle aja. Siin on aga veel üks aga! Erakonna igakuised toetusnumbrid on ise sammuti üks statistiline hinnang (enamasti saadud umbes 1000se valimiga küsitlusest)! Milline on erakondade tõeline toetus populatsioonis on kindlaks tehtav vaid kogupopulatsiooni küsitledes ja seda me teha ei saa. Sellest tulenevalt on meie küsitlustulemustes alati mingi viga sees. Mida saame teha, on hinnata toetusprotsendi vahemikku, kus tõeline toetusprotsent tõenäoliselt asetseks.
Teoreetiliselt tähendaks see olukorda, kus võtaksime populatsioonist, kus on tegelikud teotused erakondadele, lõpmatu arv juhuvalimeid ja leiaks nende valimite toetusprotsendid. Enamasti oleks tulemus sarnane tegeliku toetusega (mida me ei tea), samas teatud olukordades oleks tulemus ühele või teisele poole väga mööda. Kui me teeks nendelt juhuvalimitelt toetusprotsentide graafi, tuleks see normaaljaotusega graaf: 95% juhtumitest langeks +/-1.96 z-skoori vahele. Sellest tulenevalt saaksime hinnata, millisesse vahemikku jääb tõeline tulemus (95% tõenäosusega).
Praktikas seda teha ei saa, sest küsitlused on kallid. Küll aga on võimalik tuletada usaldusvahemik kasutades selleks standardviga. Standardviga kirjeldab, kui lähedal on valimi keskmine populatsiooni keskmisele, kasutades selleks üldjuhul variatsiooni ja valimi suurust. Standardvea leidmiseks on erinevaid meetodeid, sõltuvalt andmete iseloomust. Antud juhul on tegu proportsioonidega (erakonna toetust väljendatakse protsendina populatsioonist), mispuhul kasutame proportsiooni standardviga:
\(SE_p = \sqrt{\frac{p(1-p)}{N}}\)
kus:
Usaldusvahemiku leidmiseks tuleb korrutada standardviga soovitud z-skoori väärtusega (nt 95% usaldusvahemiku leidmiseks on see 1.96) ning liita või lahutada erakonna toetusprotsendist:
\(p \pm ci_{95} = p \pm 1.96 \cdot SE\)
Kuivõrd tegu on küllalt spetsiifilise küsimusega, siis on otstarbekas teha oma funktsioon.
ci = function(p,n=1000,z=1.96){
list(u = p + (z*(sqrt(p*(1-p)/n))),
l = p - (z*(sqrt(p*(1-p)/n))))
}
df %>%
mutate(ylemine = ci(support)[["u"]],
alumine = ci(support)[["l"]]) %>%
melt(c("party", "date"), c("support", "ylemine", "alumine"),"Toetus") %>%
ggplot() +
geom_line(aes(as.Date(date), value, color=Toetus)) +
facet_wrap(~party) +
scale_y_continuous(labels=scales::percent) +
labs(x="",y="protsent") +
theme_bw()
Lihtsamaks lugemiseks võime jooned üle siluda.
df %>%
mutate(ylemine = ci(support)[["u"]],
alumine = ci(support)[["l"]]) %>%
melt(c("party", "date"), c("support", "ylemine", "alumine"),"Toetus") %>%
ggplot() +
geom_smooth(aes(as.Date(date), value, color=Toetus), se=F, method="lm", formula= y~poly(x, 7), size=.7) +
facet_wrap(~party) +
scale_y_continuous(labels=scales::percent) +
labs(x="",y="protsent") +
theme_bw()
Alternatiivne lähenemine eelnevalt kirjeldatud meetodile oleks simuleerida erakondade toetusprotsente bootstrap meetodil. Bootstrapi idee seisneb olemasolevast valimist juhuslike valimite võtmises ning nende põhjal huvipakkuvate statistikute arvutamises (eeskätt variatsioon ja standardviga).
Kuna antud juhul on tegu proportsioonidega ning me soovime igat üksikut vaatlust eraldi simuleerida, siis võime seda teha mündiviske loogikal. Mündiviske puhul simuleerib R binaarset jaotust, mille keskmine on meie poolt antud ning jaotus on binominaalne. Korrates mündiviske katset piisava arvu kordi (antud juhul katsetame 500ga), saame arvutada simulatsiooni standardhälve ja standardvea.
Mündiviske simulatsiooniks kasutame paketti mosaic ning funktsioone do(...), rFlip(...) ja confint(...). Käesolevas katses: do(1000) kordab katset 1000 korda ning rFlip(party.n,x) viskab ühes katses münti nii mitu korda, kui palju on erakonnaga seotud vaatlusi, tõenäosusega x. confint(...) arvutab simulatsiooni usaldusvahemikud.
toetus.sim = function(x,party){
require(mosaic)
if(is.na(x)) return(NA)
party.n = nrow(df[df$party==party & !is.na(df$support),])
tmp = do(1000) * rflip(party.n,x)
ci = confint(tmp$prop)
list(u = ci[["97.5%"]],
l = ci[["2.5%"]])
}
df %>%
rowwise() %>%
mutate(ylemine = toetus.sim(support,party)[["u"]],
alumine = toetus.sim(support,party)[["l"]]) %>%
melt(c("party", "date"), c("support", "ylemine", "alumine"),"Toetus") %>%
ggplot() +
geom_smooth(aes(as.Date(date), value, color=Toetus), se=F, method="lm", formula= y~poly(x, 7), size=.7) +
facet_wrap(~party) +
scale_y_continuous(labels=scales::percent) +
labs(x="",y="protsent") +
theme_bw()
library(tidyverse)
library(effsize)
df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
paneel = read.csv("https://martenveskimae.github.io/stat-mod/ep2014_paneel.csv")
Vaatame Eesti poliitika sukka ja saabast ehk Reformierakonda ja Keskerakonda eelmise valimistsükli ajal.
df %>%
filter(party %in% c("ref", "kesk"),
year>=2011 & year<2015) %>%
ggplot() +
geom_point(aes(as.Date(date),support)) +
facet_wrap(~party) +
scale_y_continuous(labels=scales::percent) +
labs(x="") +
theme_bw()
Nende keskmine toetus on antud perioodil:
df %>%
filter(party %in% c("ref", "kesk"),
year>=2011 & year<2015) %>%
group_by(party) %>%
summarise(keskmine = mean(support,na.rm=T))
# A tibble: 2 x 2
party keskmine
<fctr> <dbl>
1 kesk 0.2444681
2 ref 0.2772340
Tekib küsimus, kas Keskerakonna keskmine toetus üldse erineb oluliselt Reformi keskmisest toetusest? Kahe keskmise võrdlemiseks saab kasutada t-testi meetodit (R-is t.test(...)), mis hindab, kas kaks keskmist pärinevad samast populatsioonist või on tegu erinevate gruppidega. Lihtsustatult öeldes, võrreldakse t-testis keskmiste erinevust kogu variatsiooniga, ning hinnatakse saadud tulemuse tõenäosust (kas tegu oli tõenäolise või ebatõenäolise tulemusega, ehk statistiliselt oluline või mitte). T-testi eelduseks on, et mõlemad valimid on normaaljaotuslikud ning sarnase variatiivsusega (viimast eeldust on võimalik rikkuda kasutades “robustset” t-testi, kus vabadusastmeid arvutatakse teisel meetodil!; R kasutab vaikimisi robustset versiooni - argument var.equal = FALSE).
T-statistiku arvutamiseks ühe valimi puhul on valemiks:
\(t=\frac{\bar{x}-\mu_o}{\frac{s}{\sqrt{N}}}\)
kus
..kahe sõltumatu valimi puhul:
\(t=\frac{\bar{x_1}-\bar{x_2}}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}\)
..ja kahe sõltuva valimi puhul:
\(t=\frac{\frac{\sum{d}}{N}}{\sqrt{\frac{\sum{d^2}-\frac{(\sum{d^2})}{N}}{N(N-1)}}}\)
kus
Teeme meetodi illustreerimiseks simulatsiooni. Ütleme, et meil on kaks sõltumatut normaaljaotusega vektorit, mõlemad 50 vaatlusega. Arvutame vektorite t-testi tulemused ja kordame katset 1000 korda (iga kord uute andmetega).
t.statistikud = replicate(1000, t.test(rnorm(50), rnorm(50))$statistic)
Vaatame saadud t-statistikute jaotust visuaalselt
ggplot() +
geom_density(aes(t.statistikud)) +
theme_bw()
Kui meie eeldused peavad paika, siis peaks tulemus lähenema t-jaotusele vabadusastmel 98 (50+50-2). Genereerime t-jaotuse antud vabadusastmetega ja arvutame 95% usaldusvahemikud. T-jaotuse leidmiseks peaksime leidma t-statistikute umbkaudsed otspunktid, näiteks käsklusega range(...), seejärel looma sümmeetrilise (sest t-jaotus on sümmeetriline!) vektori käsklusega seq(...) ning arvutama saadud vektori t-jaotusliku tiheduse etteantud vabadusastmete puhul (praegu 98), kasutades käsklust dt(...). Nendest saavad teoreetilise jaotuse x ja y atribuudid. Usaldusvahemike leidmiseks kasutame käsklust qt(..), mis arvutab teoreetilise jaotuse kvantiilid (praegu 0.025 ja 0.975).
Lisame tulemused eelnenud graafile.
r = mean(abs(range(t.statistikud)))
x = seq(-r,r,length=100)
y = dt(x,df=98)
a95 = qt(c(.025, .975), df=98)
ggplot() +
geom_density(aes(t.statistikud,color="simuleeritud t-statistikud")) +
geom_line(aes(x,y,color="teoreetiline t-jaotus")) +
geom_vline(xintercept = a95[1]) +
geom_vline(xintercept = a95[2]) +
theme_bw()
Joonte vahele jäänud tulemused olid statistiliselt ebaolulised ning servadesse sattunud tulemused statistiliselt olulised 95% usaldusnivoo puhul. Võime arvutada, mitu protsenti osutus olulisteks.
length(t.statistikud[t.statistikud<=a95[1] | t.statistikud>=a95[2]]) / length(t.statistikud)
[1] 0.043
Hindame seda esmalt ühe valimi t testiga. Leiame Reformierakonna keskmise toetuse perioodil 2011-2014 ja võrdleme seda keskerakonna toetusega.
ref.keskmine = df %>%
filter(party == "ref",
year>=2011 & year<2015) %>%
pull(support) %>%
na.omit() %>%
mean()
kesk = df %>%
filter(party == "kesk",
year>=2011 & year<2015) %>%
pull(support) %>%
na.omit()
t.test(kesk, mu=ref.keskmine)
One Sample t-test
data: kesk
t = -8.034, df = 46, p-value = 2.635e-10
alternative hypothesis: true mean is not equal to 0.277234
95 percent confidence interval:
0.2362586 0.2526775
sample estimates:
mean of x
0.2444681
Kas valitsuse ja opositsiooni toetus erineb üldse oluliselt erinevatel ajahetkedel valimistsüklis? Selleks peaksime tegema kahe sõltumatu valimi t-testi. Testime seda esialgu vaid ühel aastal, näiteks 2014.
df %>%
filter(year==2014) %>%
select(support,gov) %>%
t.test(support~gov,data=.)
Welch Two Sample t-test
data: support by gov
t = -2.3091, df = 45.844, p-value = 0.0255
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.055047311 -0.003770081
sample estimates:
mean in group 0 mean in group 1
0.2173913 0.2468000
Vaatame, kuidas t-testi p-väärtus varieerub erinevatel aastatel.
plyr::ldply(unique(df$year), function(x){
data.frame(p.value = df %>%
filter(year==x) %>%
select(support,gov) %>%
t.test(support~gov,data=.) %>%
.$p.value,
year = x)
}) %>%
ggplot() +
geom_point(aes(year,p.value)) +
geom_hline(yintercept = 0.05) +
scale_x_continuous(breaks=2007:2017) +
theme_bw()
Paneme lisaks kõrvale keskmised toetusnumbrid antud aastatel.
df %>%
group_by(year,gov=as.factor(gov)) %>%
summarise(support = mean(support,na.rm=T)) %>%
ggplot() +
geom_point(aes(year,support,color=gov)) +
geom_line(aes(year,support,color=gov)) +
scale_x_continuous(breaks=2007:2017) +
scale_y_continuous(labels=scales::percent) +
theme_bw()
Kahe sõltuva valimi t-testi erakonnatoetuse andmestikuga teha ei saa, sest meie uuringu disain ei vasta sellele. Selleks avame ühe teise andmestiku, milleks on paneeluuring enne ja pärast 2014.a. EP valimisi.
summary(paneel)
t1.date t1.evalimiste.usaldus t2.date
2014-05-20:115 Min. : 0.00 2014-06-05:222
2014-05-12:110 1st Qu.: 5.00 2014-06-03:183
2014-05-22:109 Median : 8.00 2014-06-09:173
2014-05-13:106 Mean :18.41 2014-06-04:154
2014-05-21:106 3rd Qu.:10.00 2014-06-06:131
2014-05-19:105 Max. :98.00 (Other) :139
(Other) :710 NA's :359
t2.osalusviis t2.evalimiste.usaldus
E-hääletas :141 Min. : 0.00
Ei valinud :243 1st Qu.: 5.00
Hääletas paberil:317 Median : 8.00
NA's :660 Mean :16.48
3rd Qu.:10.00
Max. :97.00
NA's :359
Vaatame näiteks, kas kurikuulus „Haldermanni case“ avaldas mõju e-valimiste usaldusele (vt näit Eksperdid: Eesti e-valimised on nii ebaturvalised, et tuleks kohe ära lõpetada). Antud asja kajastus levis laiemalt 13. mail selleks ajaks oli paneeluuringu esimeses laines juba piisavalt palju inimesi intervjueeritud, vaatame kas ja kuidas nende usaldus nägi välja enne ja kuidas pärast.
Esimene intervjueerimislaine kestis ka pärast antud juhtumit, mistõttu tuleks sealt välja filtreerida vaatlused pärast 13. maid. Ühtlasi puhastame andmed keeldunud ja puuduvatest väärtustest (kodeeritud 97 ja 98).
enne.parast = paneel %>%
mutate(t1.date = as.Date(t1.date),
usaldus.enne = case_when(
t1.date < as.Date("2014-05-13") & !(t1.evalimiste.usaldus %in% c(97:98)) ~ t1.evalimiste.usaldus),
usaldus.parast = case_when(
!(t2.evalimiste.usaldus %in% c(97:98)) ~ t2.evalimiste.usaldus)
) %>%
select(usaldus.enne, usaldus.parast)
summary(enne.parast[c("usaldus.enne", "usaldus.parast")])
usaldus.enne usaldus.parast
Min. : 0.000 Min. : 0.000
1st Qu.: 5.000 1st Qu.: 5.000
Median : 8.000 Median : 8.000
Mean : 6.767 Mean : 6.651
3rd Qu.:10.000 3rd Qu.: 9.000
Max. :10.000 Max. :10.000
NA's :1022 NA's :468
Kuidas muuta andmeid R-is juhul kui tingimusel. Kõige lihtsam variant on kasutada ifelse(...) käsklust, millesse läheb kolm argumenti:
Näiteks võtame vektori, mis sisaldab erinevaid taimi ja kodeerime selle ümber puu- ja köögiviljadeks:
taimed = c("apelsin", "kurk", "õun", "porgand", "peet")
ifelse(taimed %in% c("apelsin", "õun"), "puuvili", "köögivili")
[1] "puuvili" "köögivili" "puuvili" "köögivili" "köögivili"
Mida aga teha juhul, kui gruppe on üle kahe? If else puhul tuleks lisada iga järgnev kriteerium väära väärtussesse:
taimed = c("sibul", "till", "küüslauk", "kapsas", "piparmünt")
ifelse(taimed %in% c("sibul", "küüslauk"), "sibulköögivili",
ifelse(taimed %in% c("till", "piparmünt"), "maitseköögivili", "lehtköögivili"))
[1] "sibulköögivili" "maitseköögivili" "sibulköögivili" "lehtköögivili"
[5] "maitseköögivili"
Gruppide kasvades pole if else enam seega sobilik vahend. Paindlikum on kasutada paketist dplyr tulenevat funktsiooni case_when(), millesse tuleb sisestada vaid kaks argumenti:
Argumendid tuleb eraldada sümboliga ~ ja erinevad argumendi read komadega. Väärad tulemused märgitakse vaikimise puuduvateks väärtusteks, kuid ka neile võib anda väärtuse, mispuhul oleks tingimuseks lihtsalt TRUE:
case_when(taimed %in% c("sibul", "küüslauk") ~ "sibulköögivili",
taimed %in% c("till", "piparmünt") ~ "maitseköögivili",
TRUE ~ "lehtköögivili")
[1] "sibulköögivili" "maitseköögivili" "sibulköögivili" "lehtköögivili"
[5] "maitseköögivili"
case_when argumendi ridu saab esitada lõpmatu arv. Rohkem infot: case_when.
Nüüd saame teha kahe sõltuva valimi t-testi.
t.test(enne.parast$usaldus.enne, enne.parast$usaldus.parast, paired=T)
Paired t-test
data: enne.parast$usaldus.enne and enne.parast$usaldus.parast
t = 0.62962, df = 250, p-value = 0.5295
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.1865259 0.3618247
sample estimates:
mean of the differences
0.0876494
Eelmisest ülesandest oli näha, et statistiliselt olulist seost polnud. Kuid et tulemuses päris kindel olla, tuleks võtta arvesse veel teisigi võimalusi:
Kontrollime neid variante! Esmalt vaatame järgi, kas efekt oli lühiajaline. Üks viis selleks on kollabeerida andmed intervjueerimise päevade kaupa ja võrrelda muutust ajas.
paneel %>%
mutate(usaldus = case_when(!(t1.evalimiste.usaldus %in% 97:98) ~ t1.evalimiste.usaldus)) %>%
group_by(t1.date) %>%
summarise(keskmine = mean(usaldus, na.rm=T),
sd = sd(usaldus, na.rm=T)) %>%
ggplot() +
geom_line(aes(as.Date(t1.date), keskmine, color="Keskmine")) +
geom_line(aes(as.Date(t1.date), sd, color="Standardhälve")) +
scale_x_date(date_breaks="2 days") +
labs(x="Kuupäev", y="Usaldus (0-10)") +
theme_bw() +
theme(axis.text.x = element_text(angle=30,vjust=0.6))
Tõepoolest on näha, et usaldus langeb 13. ja 14. mai, kuid seejärel taastub.
Teiseks uurime, kas erinevus võib ilmneda mingite kolmandate tunnuste alusel. Andmestikus on meil selleks tunnus osalusviisi kohta.
summary(paneel$t2.osalusviis)
E-hääletas Ei valinud Hääletas paberil NA's
141 243 317 660
Seega tuleks teha juba varasem enne-pärast võrdlus antud gruppide lõikes. Selleks kasutame ANOVA-t, R-is aov(...):
anova.data = paneel %>%
mutate(t1.date = as.Date(t1.date),
usaldus.enne = case_when(
t1.date < as.Date("2014-05-13") & !(t1.evalimiste.usaldus %in% c(97:98)) ~ t1.evalimiste.usaldus),
usaldus.parast = case_when(
!(t2.evalimiste.usaldus %in% c(97:98)) ~ t2.evalimiste.usaldus)
) %>%
select(usaldus.enne, usaldus.parast, t2.osalusviis)
aov(usaldus.enne-usaldus.parast ~ t2.osalusviis, data=anova.data) %>% summary()
Df Sum Sq Mean Sq F value Pr(>F)
t2.osalusviis 2 35.2 17.624 3.743 0.0257 *
Residuals 170 800.4 4.708
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1188 observations deleted due to missingness
Ka siin selgub, et esineb statistiliselt oluline seos. Paraku ei anna anova otsest aimdust selle kohta, milliste gruppide vahel seos on. Hea on seega vaadata jaotusi graafiliselt, näiteks karpdiagrammi abil.
anova.data %>%
na.omit() %>%
ggplot() +
geom_boxplot(aes(t2.osalusviis, usaldus.enne-usaldus.parast)) +
theme_bw()
Pilt viitab justkui sellele, et paberhääletajatel jaotusid erinevused võrdselt ümber nulli, samas e-hääletajate usaldus oli langenud ja mittehääletajatel tõusnud.
Tihti huvitab meid statistikas efekti suurus. St, kui suurt mõju üks tunnus teisele avaldab. Keskmiste võrdluse puhul loetakse efekti suuruseks erinevust erinevust kahe keskmise vahel, jagatuna andmete standardhälvega:
\(d=\frac{\bar{x_1}-\bar{x_2}}{s}\)
Antud efekti suurust nimetatakse Coheni d-ks. Coheni d leidmiseks kasutame paketti effsize ja funktsiooni cohen.d(...), mis jäljendab oma loogikalt t-testi t.test(...). (Antud funktsioon kasutab siiski eelnevalt toodud valemi robustsemat edasiarendust, kuid selle loogika on siiski sama!)
df %>%
filter(party %in% c("kesk", "ref")) %>%
cohen.d(support ~ party, data=.)
Cohen's d
d estimate: 0.536684 (medium)
95 percent confidence interval:
inf sup
NaN NaN
library(tidyverse)
df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
paneel = read.csv("https://martenveskimae.github.io/stat-mod/ep2014_paneel.csv")
Tulpdiagramm on sobilik juhul, kui soovime võrrelda kategoorilise tunnuse gruppe ühe või väikese arvu tunnuste lõikes. Juhul, kui tunnuseid on mitu, võib tunnuste väärtused asetada üksteise otsa või üksteise kõrvale (lähemalt all). Antud diagrammi eelis seisneb võimaluses hinnata lihtsa vaevaga tulpade erinevusi tulpade kõrgustes (sellest tulenevalt tuleks võrreldavate tunnuste arv hoida pigem madalal).
Püüame näiteks graafida erakondade keskmised toetusnumbrid aastate lõikes. Tulpdiagrammi loomiseks on funktsioon geom_col(...). Teeme tulpdiagrammi, kus x-teljel on aastad, y-teljel iga erakondade keskmine toetus ning tulba värv on seotud erakonnaga.
df %>%
group_by(party, year) %>%
summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
ggplot() +
geom_col(aes(year, keskmine.toetus, fill=party))
Antud juhul käsitletakse x-teljel olevat tunnust year arvtunnusena, mistõttu on ggplot ise otsustanud, millise vahega väärtusi joonisel kuvada. Praegu võiks kuvada kõik aastad. Kõige lihtsam variant on kasutada käsklust as.factor(year), mis muudab tunnuse kategooriliseks (võib kasutada nii ggploti sees, kui ka enne seda).
df %>%
group_by(party, year = as.factor(year)) %>%
summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
ggplot() +
geom_col(aes(year, keskmine.toetus, fill=party))
Joonisel on keeruline võrrelda tulpade keskel olevaid erakondi. Parem oleks asetada erinevad erakonnad üksteise kõrvale. Selleks kasutame argumenti position: geom_col(position="dodge").
df %>%
group_by(party, year = as.factor(year)) %>%
summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
ggplot() +
geom_col(aes(year, keskmine.toetus, fill=party), position="dodge")
Parem! Järgmiseks oleks hea muuta joonise värve üksteisest selgemalt eristuvaks. Muudame ära ka legendi pealkirja. Selleks peaksime muutma fill skaala argumente scale_fill_brewer(palette="Accent", name="Erakonnad") (värvide käsitsi sisestamiseks tuleks kasutada scale_fill_manual(values=c(...))).
R kasutab värvide genereerimiseks Color Brewerit (pakett RColorBrewer). Rohkem infot: R Color Cheatsheet
df %>%
group_by(party, year = as.factor(year)) %>%
summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
ggplot() +
geom_col(aes(year, keskmine.toetus, fill=party), position="dodge") +
scale_fill_brewer(palette="Accent", name="Erakonnad")
Veel tuleks muuta y-telje skaala protsentideks, kasutades selleks paketti scales scale_y_continuous(labels=scales::percent). Muudame ära ka telgede nimetused ning lisame pealkirja ja viite käsklusega labs(x="Aasta", y="Toetus",title="Erakondade keskmine toetus",caption="Kantar Emor küsitlusandmetel").
df %>%
group_by(party, year = as.factor(year)) %>%
summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
ggplot() +
geom_col(aes(year, keskmine.toetus, fill=party), position="dodge") +
scale_fill_brewer(palette="Accent", name="Erakonnad") +
scale_y_continuous(labels=scales::percent) +
labs(x="Aasta", y="Toetus",title="Erakondade keskmine toetus",caption="Kantar Emor küsitlusandmetel")
Viimaks muudame joonise tausta valgeks juba tuttava käsklusega theme_bw().
df %>%
group_by(party, year = as.factor(year)) %>%
summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
ggplot() +
geom_col(aes(year, keskmine.toetus, fill=party), position="dodge") +
scale_fill_brewer(palette="Accent", name="Erakonnad") +
scale_y_continuous(labels=scales::percent) +
labs(x="Aasta", y="Toetus",title="Erakondade keskmine toetus",caption="Kantar Emor küsitlusandmetel") +
theme_bw()
Histogramm on sobilik arvtunnuse jaotuse visualiseerimiseks. Histogrammi kuvamiseks on ggplotis käsklus geom_histogramm(...). Histogrammi y-telg on vaikimisi sagedus, ent selle võib kirjutada soovi korral üle. Histogrammil on võimalik seadistada vahemike laiust või hulka (argumendid bin ja binwidth).
df %>%
ggplot() +
geom_histogram(aes(support), binwidth=0.02) +
scale_x_continuous(labels=scales::percent) +
labs(x="Toetus", y="Sagedus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
theme_bw()
Iga erakonna kohta eraldi joonise kuvamiseks võime kasutada käsklusi facet_wrap(~party) või facet_grid(~party).
df %>%
ggplot() +
geom_histogram(aes(support), binwidth=0.02) +
facet_wrap(~party) +
scale_x_continuous(labels=scales::percent) +
labs(x="Toetus", y="Sagedus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
theme_bw()
Arvtunnuse jaotuse näitamiseks on võimalik kasutada lisaks histogrammile karpdiagrammi geom_boxplot(...) ja viiuldiagrammi geom_violin(...). Karpdiagrammile võib hajutatult kuvada ka üksikud vaatlused käsklusega geom_jitter(...) (mõistlik on muuta kuvatud vaatlused läbipaistvaks argumendiga alpha ning koondada need kitsamale alale argumendiga width).
df %>%
ggplot() +
geom_boxplot(aes(party,support)) +
geom_jitter(aes(party,support), alpha=0.2, width=0.2) +
scale_y_continuous(labels=scales::percent) +
labs(x="Erakond", y="Toetus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
theme_bw()
Viiuldiagramm annab veelgi täpsema ülevaate jaotusest. Viiuldiagrammile võib lisada kvantiilide jooned argumendiga draw_quantiles.
df %>%
ggplot() +
geom_violin(aes(party,support), draw_quantiles=c(.25,.5,.75)) +
scale_y_continuous(labels=scales::percent) +
labs(x="Erakond", y="Toetus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
theme_bw()
Hajuvusdiagrammi kasutatakse kahe arvtunnuse omavahelise seose visualiseerimiseks.
df %>%
ggplot() +
geom_point(aes(support,gdpgr_yearly)) +
scale_x_continuous(labels=scales::percent) +
facet_wrap(~party) +
labs(x="Toetus", y="SKT aastane muutus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
theme_bw()
Joondiagramm, nagu hajuvusdiagrammgi, on sobilik kahe arvtunnuse omavahelise seose visualiseerimiseks. Joondiagrammi kasutatakse rohkem siis, kui üks tunnustest on seotud ajaga.
df %>%
ggplot() +
geom_line(aes(as.Date(date),support)) +
scale_y_continuous(labels=scales::percent) +
facet_wrap(~party) +
labs(x="", y="Toetus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
theme_bw()
Küsitlustes kasutatakse tihti likerti skaalasid, mida visualiseeritakse enamasti tulpdiagrammi abil (skaala on protsentides).
Kasutame näiteks EP 2014 andmestikust usalduse tunnust.
library(reshape2)
usaldus = paneel %>%
mutate(enne = case_when(!(t1.evalimiste.usaldus %in% 97:98) ~ t1.evalimiste.usaldus),
parast = case_when(!(t2.evalimiste.usaldus %in% 97:98) ~ t2.evalimiste.usaldus)) %>%
select(enne,parast) %>%
melt() %>%
na.omit() %>%
group_by(variable,value = as.factor(value)) %>%
summarise(count = n()) %>%
mutate(pc = count/sum(count[variable==variable]))
usaldus
# A tibble: 22 x 4
# Groups: variable [2]
variable value count pc
<fctr> <fctr> <int> <dbl>
1 enne 0 185 0.15651438
2 enne 1 17 0.01438240
3 enne 2 29 0.02453469
4 enne 3 44 0.03722504
5 enne 4 32 0.02707276
6 enne 5 85 0.07191201
7 enne 6 42 0.03553299
8 enne 7 88 0.07445008
9 enne 8 185 0.15651438
10 enne 9 176 0.14890017
# ... with 12 more rows
Graafime tulemused tulpdiagrammina.
usaldus %>%
ggplot() +
geom_col(aes(variable,pc,fill=value))
Antud joonist saab muuta loetavamaks vahetades värviskaalat, pöörates joonise telgi, muutes telgede nimesid jne.
usaldus %>%
ggplot() +
geom_col(aes(variable,pc,fill=reorder(value,desc(value))),width=0.7) +
scale_fill_brewer(name="Usaldus",palette="Paired") +
scale_y_continuous(labels=scales::percent) +
scale_x_discrete(labels=c("Enne","Pärast")) +
labs(x="", y="Osakaal", title="E-valimiste usaldus",caption="EP 2014 valimiste küsitlusandmed") +
coord_flip() +
theme_minimal()
Alternatiiv on kasutada spetsiaalset likerti skaalade visualiseerimiseks mõeldud funktsiooni likert(...) paketist HH.
usaldus.2 = usaldus %>%
dcast(variable ~ value, value.var="count")
rownames(usaldus.2) = c("Enne", "Pärast")
HH::likert(usaldus.2,
as.percent=T,
ylab.right="Vastanute arv",
xlab="Osakaal",
main="E-valimiste usaldus",
sub="")
library(tidyverse)
library(psych)
library(PerformanceAnalytics)
df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
paneel = read.csv("https://martenveskimae.github.io/stat-mod/ep2014_paneel.csv")
Vaatame lähemalt kahe tunnuse vahelist korrelatsiooni, antud juhul Pearsoni korrelatsiooni. Pearsoni korrelatsioon sobib juhul, kui mõlema tunnuse puhul on tegu arvtunnustega.
Pearsoni korrelatsioon on iseloomult lineaarne, st kahe tunnuse vahele üritatakse sobitada lineaarset seost. Matemaatiliselt on korrelatsioon kahe tunnuse kovariatsioon, skaleerituna nende standardhälbega:
\(r=\frac{\sum\limits_{i=1}^n(x_{i}-\bar{x})(y_{i}-\bar{y})}{(n-1)s_{x}s_{y}}\)
Et anda hinnang korrelatsioonikordaja statistilisele olulisusele, võime arvutada korrelatsiooni t-statistiku väärtuse ja võrrelda seda teoreetilise t-jaotusega.
\(t=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}\)
Arvutame käsitsi kahe juhusliku vektori korrelatsiooni (eeldatavalt on korrelatsioon juhuslikkuse tõttu 0). Esmalt tekitame kaks vektorit ning vaatame võimalikku visuaalselt.
n = 100
x = rnorm(n)
y = rnorm(n)
plot(x,y)
Kui kahe vektori vahel on tõepoolest seos, siis on see juhuslik. Hindame seost korrelatsiooni abil.
xbar = mean(x)
ybar = mean(y)
sx = sd(x)
sy = sd(y)
r = sum((x-xbar)*(y-ybar)) / ((n-1)*sx*sy)
r
[1] 0.1136141
Kas tulemus on statistiliselt oluline?
t = (r*sqrt(n-2)) / sqrt(1-r^2) # t-väärtus
p = 2*pt(-abs(t),df=199) # p-väärtus
# Olulisuse visualiseerimine (sama, mis t-testi puhul)
tdistx = seq(-3,3,length=300)
tdisty = dt(tdistx,df=199)
a95 = qt(c(.025, .975), df=199)
ggplot() +
geom_line(aes(tdistx,tdisty)) +
geom_vline(xintercept = a95[1]) +
geom_vline(xintercept = a95[2]) +
geom_vline(xintercept = t, color="red") +
theme_bw()
Arvutuste tulemusel on:
Õnneks teeb kogu töö ära funktsioon corr.test(...) paketist psych:
corr.test(cbind(x,y))
Call:corr.test(x = cbind(x, y))
Correlation matrix
x y
x 1.00 0.11
y 0.11 1.00
Sample Size
[1] 100
Probability values (Entries above the diagonal are adjusted for multiple tests.)
x y
x 0.00 0.26
y 0.26 0.00
To see confidence intervals of the correlations, print with the short=FALSE option
Tulemus on sama!
Majandusliku hääletamise teooriast teame, et hääletamine on tihti seotud majanduse käekäiguga - kui majandusel läheb hästi, on koalitsiooni toetus kõrge, kui halvasti, on opositsiooni toetus kõrge. Vaatame näiteks, kuidas on seotud Reformierakonna toetus tööpuuduse, SKT kasvu ja inflatsiooniga.
ref = df %>%
filter(party == "ref") %>%
select(support, unemp, gdpgr_yearly, inflats_m) %>%
na.omit()
cor(ref)
support unemp gdpgr_yearly inflats_m
support 1.0000000 -0.0712523 0.3312793 0.3958776
unemp -0.0712523 1.0000000 -0.1915624 -0.2088902
gdpgr_yearly 0.3312793 -0.1915624 1.0000000 0.3311977
inflats_m 0.3958776 -0.2088902 0.3311977 1.0000000
Korrelatsioonide statistilise olulisuse jaoks tuleks kasutada corr.test(...) funktsiooni paketis psych.
corr.test(ref)
Call:corr.test(x = ref)
Correlation matrix
support unemp gdpgr_yearly inflats_m
support 1.00 -0.07 0.33 0.40
unemp -0.07 1.00 -0.19 -0.21
gdpgr_yearly 0.33 -0.19 1.00 0.33
inflats_m 0.40 -0.21 0.33 1.00
Sample Size
[1] 93
Probability values (Entries above the diagonal are adjusted for multiple tests.)
support unemp gdpgr_yearly inflats_m
support 0.0 0.50 0.01 0.00
unemp 0.5 0.00 0.13 0.13
gdpgr_yearly 0.0 0.07 0.00 0.01
inflats_m 0.0 0.04 0.00 0.00
To see confidence intervals of the correlations, print with the short=FALSE option
Vaikimisi annab antud funktsioon nii korrelatsioonimaatriksi kui p-väärtused. Kui soovida vaid ühte, tuleks käskluse lõppu lisada kas $r või $p korrelatsioonide või p-väärtuste saamiseks.
corr.test(ref)$r
support unemp gdpgr_yearly inflats_m
support 1.0000000 -0.0712523 0.3312793 0.3958776
unemp -0.0712523 1.0000000 -0.1915624 -0.2088902
gdpgr_yearly 0.3312793 -0.1915624 1.0000000 0.3311977
inflats_m 0.3958776 -0.2088902 0.3311977 1.0000000
Alati on mõistlik hinnata pilti graafiliselt. Korrelatsiooni puhul kasutatakse tihti spetsiaalseid korrelatsiooni visualiseeringuid, mida on keeruline ggplotis luua. Küll aga on mitmeid pakette, mis antud probleemi lahendavad. Näiteks võib kasutada funktsiooni chart.Correlation(...) paketist PerformanceAnalytics.
chart.Correlation(ref)
Juhul, kui tunnused on järjestikused, on mõistlik kasutada Spearmani korrelatsiooni, mis ei kasuta korrelatsiooni arvutamiseks algtunnused, vaid annab neile astmeväärtuse (rank). Spearmani korrelatsiooni läbiviimiseks võib kasutada juba tuttavat EP 2014 valimiste andmestikku ning usalduse tunnuseid.
usaldus = paneel %>%
mutate(enne = case_when(!(t1.evalimiste.usaldus %in% 97:98) ~ t1.evalimiste.usaldus),
parast = case_when(!(t2.evalimiste.usaldus %in% 97:98) ~ t2.evalimiste.usaldus)) %>%
select(enne,parast) %>%
na.omit()
corr.test(usaldus, method="spearman")
Call:corr.test(x = usaldus, method = "spearman")
Correlation matrix
enne parast
enne 1.00 0.75
parast 0.75 1.00
Sample Size
[1] 835
Probability values (Entries above the diagonal are adjusted for multiple tests.)
enne parast
enne 0 0
parast 0 0
To see confidence intervals of the correlations, print with the short=FALSE option
chart.Correlation(usaldus, method="spearman")
library(tidyverse)
library(reshape2)
df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
Mitmed mudelid R-is, sh lineaarse regressiooni mudel, küsivad sisendit valemi kujul. Meenutame näiteks lineaarse regressiooni valemit:
\(\hat{y} = \alpha + \beta X + \epsilon\)
kus:
Lineaarse regressiooni jooksutamiseks oleks vaja sisendiks anda nii y kui ka x tunnused, misjärel arvutatakse välja \(\alpha\) ja \(\beta\) väärtused.
Jätkame siinkohal korrelatsioonanalüüsis uuritud tunnustega - erakonna toetus ja töötus. Näitena kasutame taas Reformierakonda. Valemi kujul näeks meie mudel välja järgnev:
\(toetus_{ref} = \beta_0 + \beta_1 töötus + \epsilon\)
Mudeli spetsifitseerimisel näeks see välja toetus ~ töötus, kus ~ sümbol eristab sõltuvat ja sõltumatuid tunnuseid. Lineaarse regressiooni jooksutamiseks kasutame funktsiooni lm(...), kuhu sisendiks läheb valem ning andmestiku nimetus (data=...).
Filtreerime andmestikust Reformierakonnaga seotud vaatlused ning kasutame support ja unemp tunnuseid. Enne analüüsini jõudmist viime veel support tunnuse 0 ja 100 vahele, et ühtlustada erinevate tunnuste skaalad ning lihtsustada hilisemat tulemuste tõlgendamist. Lisaks arvutame aastase inflatsiooni, mida kasutame mitmese lineaarse regressiooni näites.
df = df %>%
mutate(support=support*100,
inflats_a =((inflats_m_basevalue / inflats_prevyear_base) - 1)*100)
ref = df %>% filter(party=="ref")
ref_mudel = lm(support ~ unemp, data=ref)
ref_mudel
Call:
lm(formula = support ~ unemp, data = ref)
Coefficients:
(Intercept) unemp
27.3345 0.2048
Mudeli jooksutamisel kuvatakse \(\beta_0\) (Intercept) ja \(\beta_1\) (unemp), antud juhul väärtustega 27.335 ja 0.205. Kuna tunnused on nüüd samal protsentuaalsel skaalal, siis saame järeldada, et töötuse 1% suurenedes kasvab Reformierakonna toetus 0.2%. Arvestades, et Reformierakond on olnud pikalt valitsuses, siis võiks eeldada hoopis negatiivset seost - töötuse kasvades langeb valitseva erakonna toetus. Mudeliga tuleks seega veel tööd teha, ent esmalt uurime teisi mudeliga seotud statistikuid.
summary(ref_mudel)
Call:
lm(formula = support ~ unemp, data = ref)
Residuals:
Min 1Q Median 3Q Max
-12.6659 -5.0756 0.2317 3.7196 16.6618
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 27.3345 1.6422 16.645 <2e-16 ***
unemp 0.2048 0.1615 1.268 0.208
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.772 on 107 degrees of freedom
(15 observations deleted due to missingness)
Multiple R-squared: 0.01481, Adjusted R-squared: 0.005599
F-statistic: 1.608 on 1 and 107 DF, p-value: 0.2075
summary(...) funktsioon annab meile lisaks \(\beta\) väärtuste standard vead ja p-väärtused, ning \(R^2\) väärtused. Näeme, et töötuse tunnus ei ole statistiliselt oluline (\(\beta_1\) ei erine oluliselt nullist) ning mudeli seletusvõime on vaid 1.5% (Multiple R-squared). Ühtlasi näeme viimaselt realt, et kõiki sõltumatuid tunnuseid arvesse võttes pole mudel statistiliselt oluline (p-väärtus: 0.208). Kuna antud juhul kaasasime vaid ühe sõltumatu tunnuse, ei paku see uut teadmist, ent keerulisemate mudelite puhul, kus on mitmeid sõltumatuid tunnuseid, tasub seda jälgida.
Regressioonanalüüsis hinnatakse \(\beta\) väärtuste statistilist olulisust t-testi abil, hinnates, kas väärtus on oluliselt erinev nullist. Sellest tulenevalt, mida lähemal on \(\beta\) nullile, seda suurem kipub olema tema p-väärtus. Seega, väikeste \(\beta\) väärtuste puhul tuleks pidada meeles, et statistiliselt olulise tulemuse saamise tõenäosus on madalam.
Alljärgnevalt viime läbi sumulatsiooni ning uurime \(\beta_0\) ja \(\beta_1\) varieeruvust erinevate valimi suuruste puhul. Hindame seejuures tüüp I viga ning näeme, et see vastab valitud usaldusnivoole.
Populatsiooni parameetrite suurusi võib muuta, et piltlikustada eelpool kirjeldatud nähtust.
N = 100000 # Populatsiooni suurus
n = 1000 # Valimi ülemine piir
x = rnorm(N,50,15) # X
e = rnorm(N,0,10) # Juhuslik viga
b0 = 2.5 # beta 0
b1 = 5 # beta 1
y = b0 + (b1*x) + e # Y
# Loome andmetabeli x ja y väärtustega
pop = data.frame(x = x, y = y)
# Loome funktsiooni, mis väljastab regressiooni tulemused vastavalt valimi suurusele
reg = function(data,n){
df = data[sample(1:nrow(data), n),]
m = lm(y~x,data=df)
data.frame(b0 = m$coefficients[[1]],
b0l = confint(m)[1,1],
b0u = confint(m)[1,2],
b0p = summary(m)$coefficients[,4][[1]],
b0se = summary(m)$coefficients[,2][[1]],
b1 = m$coefficients[[2]],
b1l = confint(m)[2,1],
b1u = confint(m)[2,2],
b1p = summary(m)$coefficients[,4][[2]],
b1se = summary(m)$coefficients[,2][[2]],
n=n)
}
# Jooksutame regressiooni erinevate valimi suurustega
r = lapply(10:n,function(x)reg(pop,x)) %>% do.call(rbind.data.frame,.)
# Hindame, kas saadud tulemustes esines tüüp I viga (false positive)
b0_pop = reg(pop,N)$b0
b1_pop = reg(pop,N)$b1
r = r %>%
mutate(b0_typeI = ifelse((b0l > b0_pop | b0u < b0_pop) & b0p <= .05, TRUE, FALSE),
b1_typeI = ifelse((b1l > b1_pop | b1u < b1_pop) & b1p <= .05, TRUE, FALSE))
# Graafime tulemused
r %>%
select(n,b0_typeI,b1_typeI,b0,b1) %>%
melt(c("n","b0_typeI","b1_typeI")) %>%
melt(c("n","variable","value"),variable.name="TF",value.name="type_I") %>%
ggplot() +
geom_point(aes(n,value,color=type_I),shape=1,alpha=.5) +
geom_hline(data=data.frame(variable=c("b0","b1"),
value=c(b0_pop,b1_pop)),
aes(yintercept=value)) +
scale_color_manual(values=c("steelblue","red")) +
facet_wrap(~variable) +
theme_bw()
r %>%
ggplot() +
geom_point(aes(n,b0p,color=b0_typeI),shape=1) +
geom_hline(yintercept=0.05) +
scale_color_manual(values=c("steelblue","red")) +
theme_bw()
# Arvutame, mitmel protsendil juhtudest esines tüüp I viga (eeldada võiks ligemale 5%)
sum(r$b0_typeI==TRUE)/nrow(r)
[1] 0.02825429
sum(r$b1_typeI==TRUE)/nrow(r)
[1] 0.06559031
Tihti kasutatakse mudeli hindamiseks veel ruutkeskmist hälvet (mean squared error ehk MSE), mis väljendab, kui hästi suudab mudel ennustada sõltuva tunnuse väärtuseid. MSE leidmiseks tuleb lahutada õigetest väärtustest ennustatud väärtused (ehk arvutada ennustusjäägid), võtta need ruutu ning leida saadud jada keskmine. Seda on lihtne teha, kuivõrd mudel pakub ennustusjääke $residuals abil.
mean(ref_mudel$residuals^2)
[1] 45.01441
Kui soovime näha mudeli ennustatud väärtusi, siis on selleks funktsioon predict(...), kuhu läheb argumendiks mudel ning soovi korral andmestik uute andmetega.
yhat = predict(ref_mudel)
head(yhat)
1 2 3 4 5 6
28.42012 28.33819 28.33819 28.33819 28.17432 28.17432
Sarnase mudeli võib teha ka teiste erakondade kohta. Võtame näiteks Keskerakonna.
kesk = df %>% filter(party=="kesk")
kesk_mudel = lm(support ~ unemp, data=kesk)
summary(kesk_mudel)
Call:
lm(formula = support ~ unemp, data = kesk)
Residuals:
Min 1Q Median 3Q Max
-7.1163 -2.3599 -0.3058 2.1375 9.5876
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.53591 0.82636 29.692 <2e-16 ***
unemp 0.11844 0.08128 1.457 0.148
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.408 on 107 degrees of freedom
(15 observations deleted due to missingness)
Multiple R-squared: 0.01946, Adjusted R-squared: 0.0103
F-statistic: 2.124 on 1 and 107 DF, p-value: 0.148
Tulemus on suuresti sama.
Mudeli tulemuste mugavaks esitamiseks on loodud pakett stargazer. Stargazer võimaldab esitada nii lineaarse regressiooni kui mitmete teiste mudelite tulemusi kujul, mis on aktsepteeritav enamikes teadusajakirjades. Stargazeri võimalustega tutvumiseks saab infot järgnevalt lingilt: Stargazer cheatsheet.
library(stargazer)
# type võimalikud väärtused on "text", "html" ja "latex", vastavalt soovitud väljundi formaadile (text, html ja pdf)
stargazer(ref_mudel, kesk_mudel,
title="Linear regression table",
covariate.labels="Unemployment",
dep.var.labels="Support",
column.labels=c("Ref", "Kesk"),
model.numbers=F,
no.space=T,
star.cutoffs=c(0.05,0.01,0.001),
type="html")
| Dependent variable: | ||
| Support | ||
| Ref | Kesk | |
| Unemployment | 0.205 | 0.118 |
| (0.162) | (0.081) | |
| Constant | 27.335*** | 24.536*** |
| (1.642) | (0.826) | |
| Observations | 109 | 109 |
| R2 | 0.015 | 0.019 |
| Adjusted R2 | 0.006 | 0.010 |
| Residual Std. Error (df = 107) | 6.772 | 3.408 |
| F Statistic (df = 1; 107) | 1.608 | 2.124 |
| Note: | p<0.05; p<0.01; p<0.001 | |
Lineaarset regressiooni kasutades ei pea piirduma vaid ühe tunnuse lisamisega, vaid see võimaldab hinnata mitmeid tunnuseid korraga (oluline eelis t-testi ja korrelatsiooni ees). Lisame varasemale mudelile SKT kasvu ja inflatsiooni tunnused:
\(toetus_{ref} = \beta_0 + \beta_1 töötus + \beta_2 SKT + \beta_3 inflatsioon + \epsilon\)
mis näeks R-i valemina välja support ~ unemp + gdpgr_yearly + inflats_m. Lihtsustamise mõttes piirame vaatlused eelviimase valitsustsükliga, et kontrollida koalitsiooni/opositsiooni mõju ning jooksutame uue mudeli lisatunnustega.
ref = ref %>% filter(year>=2011, year<2015)
mlr_mudel = lm(support ~ unemp + gdpgr_yearly + inflats_a, data=ref)
summary(mlr_mudel)
Call:
lm(formula = support ~ unemp + gdpgr_yearly + inflats_a, data = ref)
Residuals:
Min 1Q Median 3Q Max
-6.2222 -2.3461 0.6219 1.5866 11.2140
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 29.9480 3.8262 7.827 8.35e-10 ***
unemp -0.8965 0.5845 -1.534 0.132
gdpgr_yearly 1.8274 0.3438 5.315 3.59e-06 ***
inflats_a -0.3587 0.4936 -0.727 0.471
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.535 on 43 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.4927, Adjusted R-squared: 0.4573
F-statistic: 13.92 on 3 and 43 DF, p-value: 1.761e-06
Tulemused näitavad, et SKT kasv on statistiliselt oluline ning märkimisväärse mõjuga. Samuti kasvas mudeli seletusvõime 46%ni (Adjusted R-squared) ning mudel ise on statistiliselt oluline.
mean(mlr_mudel$residuals^2)
[1] 11.43551
Mudeli ruutkeskmine hälve on samuti oluliselt vähenenud võrreldes eelnevaga.
Kui võtta eesmärgiks saada võimalikult hea mudeli sobivus, võib kaaluda sõltumatute tunnuste transformeerimist, eeldusel, et see on teoreetiliselt põhjendatud. Seega tuleks küsida, milline funktsionaalne vorm võiks antud sõltumatutel tunnustel peale lineaarse vormi olla. Visualiseerime seosed.
ref %>%
melt("support",c("unemp","gdpgr_yearly","inflats_a")) %>%
ggplot() +
geom_point(aes(value,support)) +
coord_fixed() +
facet_wrap(~variable,scales="free",ncol=2)
Antud juhul võib jätta töötuse lineaarseks, SKT kasvu puhul katsetada logaritmimist ning inflatsiooni puhul kasutada ruutseost. Loome uued transformeeritud tunnused, ning vaatame, kas seos erakonna toetusega on muutunud lineaarsemaks.
ref = ref %>%
mutate(gdpgr_yearly_log = log(gdpgr_yearly),
inflats_a_squared = inflats_a^2)
ref %>%
melt("support",c("unemp","gdpgr_yearly_log","inflats_a_squared")) %>%
ggplot() +
geom_point(aes(value,support)) +
coord_fixed() +
facet_wrap(~variable,scales="free",ncol=2)
Tundub, et on olnud abi ning võime teha uue mudeli transformeeritud tunnustega.
mlr_mudel_tr = lm(support ~ unemp + gdpgr_yearly_log + inflats_a + inflats_a_squared, data=ref)
summary(mlr_mudel_tr)
Call:
lm(formula = support ~ unemp + gdpgr_yearly_log + inflats_a +
inflats_a_squared, data = ref)
Residuals:
Min 1Q Median 3Q Max
-6.9753 -1.9009 -0.3334 1.8826 10.5300
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 29.0877 3.8250 7.605 2.01e-09 ***
unemp -0.3867 0.5140 -0.752 0.45602
gdpgr_yearly_log 3.0627 0.9127 3.356 0.00169 **
inflats_a -2.5772 1.2790 -2.015 0.05033 .
inflats_a_squared 0.5346 0.2694 1.984 0.05379 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.481 on 42 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.5195, Adjusted R-squared: 0.4738
F-statistic: 11.35 on 4 and 42 DF, p-value: 2.461e-06
mean(mlr_mudel_tr$residuals^2)
[1] 10.83034
Mudeli parameetrites on näha väikest arengut. Edasi tuleks anda põhjalikum hinnang mudeli toimimisele.
Kõige kiirema ülevaate võimalikest probleemidest saab ennustusjääkide visualiseerimisel.
plot(mlr_mudel_tr, which=1) # Rohkem joonised ennustusjääkidest saab jättes 'which' täpsustamata
Ennustusjääkide puhul soovime näha, et need oleksid jaotunud ühtlase hajususega ümber nulli. Antud juhul on näha, et ennustusjäägid on siiski teatud määral süstemaatiliselt kaldu, mis viitab probleemidele mudeli funktsionaalses vormis. Samas heteroskedastiivsust ei paista eriti olevat.
Mudeli funktsionaalse vormi hindamiseks kasutatakse tihti Ramsey RESET-testi, mis lisab mudelile regressorite ruutseosed ning hindab, kas see on algmudelist parem. Kasutame selleks funktsiooni resettest(...) paketist lmtest.
library(lmtest)
resettest(mlr_mudel, power=2, type="fitted")
RESET test
data: mlr_mudel
RESET = 6.2981, df1 = 1, df2 = 42, p-value = 0.01603
resettest(mlr_mudel_tr, power=2, type="fitted")
RESET test
data: mlr_mudel_tr
RESET = 2.4592, df1 = 1, df2 = 41, p-value = 0.1245
Antud juhul viitavad tulemused, et ilma transformeeritud tunnusteta mudelil oli halvem funktsionaalne vorm kui transformeeritud tunnustega mudelil. Transformeeritud tunnustega mudel ei muutu oluliselt paremaks, kui lisada ruutseostega tunnused. Samas võib RESET testi tulemus olla oluline kuupseoste (argument power=2:3) korral, kuivõrd visuaalselt oli näha ennustusjääkide teatud süstemaatilist kõikumist ümber nulli.
Heteroskedastiivsuse hindamiseks on populaarsemad White test ja Breusch-Pagan test, millest võib antud juhul katsetada viimasega. Breusch-Pagani testi puhul sobitatakse uus mudel algmudeli ennustusjääkidele ning hinnatakse, kas sealt ilmnevad olulised seosed. lmtest paketist on selleks funktsioon bptest(...).
bptest(mlr_mudel_tr)
studentized Breusch-Pagan test
data: mlr_mudel_tr
BP = 5.2184, df = 4, p-value = 0.2656
Kuna antud juhul osutus test statistiliselt ebaoluliseks, võib järeldada, et heteroskedastiivsusega olulisi probleeme pole.
Viimaks võime kontrollida, ega uute tunnuste lisamisel pole mudelis ilmnenud multikollineaarsust. Multikollineaarsuse korral on sõltumatud tunnused on omavahel tugevas korrelatsioonis, mistõttu on raske eristada üksikute tunnuste mõju (korrelatsioonist tingituna, kui üks väärtus kasvab, kasvavad ka teised väärtused ning üksikud mõjud pole enam eristatavad). Sellest tulenevalt pole hinnatud parameetrid (\(\beta\) väärtused) enam usaldatavad. Kuigi teatud korrelatsioon esineb tunnuste vahel peaaegu alati, on oluline, et sõltumatute tunnuste omavaheline korrelatsioon ei ületaks sõltumatu ja sõltuva tunnuse korrelatsiooni.
Regressioonanalüüsi korral on multikollineaarsuse hindamiseks kasutusel VIF (Variance Inflation Factor) test. VIF test hindab, kui suur oli sõltumatu tunnuse \(\beta_k\) variatsioon teiste sõltumatute tunnuste kontekstis (\(Var(\beta_k)\)) ja kui suur isoleeritult (\(Var_{min}(\beta_k)\)). Seejärel jagab ta esimese teisega \(\frac{Var(\beta_k)}{Var_{min}(\beta_k)} = \frac{1}{1-R_k^2}\). Üldjuhul leitakse, et kui tulemus on üle kümne, võib mudelis esineda multikollineaarsuse probleem.
VIF testi funktsiooni vif(...) saame paketist HH. Kuna HH ja tidyiverse pakettides on osa funktsioone sama nimega, siis ei hakka praegu HH-d eraldi aktiveerima, vaid rakendame sealt ainult hetkel soovitud vif funktsiooni.
HH::vif(mlr_mudel_tr)
unemp gdpgr_yearly_log inflats_a inflats_a_squared
4.155656 2.647484 20.732871 27.740997
Antud mudeli puhul on probeelm inflatsiooni ja inflatsioon ruudus tunnustega, mis on samas ootuspärane tulemus (võimegi eeldada tugevat seost antud kahe tunnuse vahel)
HH::vif(mlr_mudel)
unemp gdpgr_yearly inflats_a
5.210994 3.358901 2.994678
Transformeerimata tunnustega mudelis sellist probleemi ei esinenud ning ilmselt oleksid selle tulemused kokkuvõttes usaldusväärsemad.
library(tidyverse)
library(stargazer)
library(ResourceSelection)
library(caret)
RK_2015 = read.csv("https://martenveskimae.github.io/stat-mod/RK_2015.csv")
Eelnevalt käsitletud lineaarse regressiooni mudeli üheks peamiseks eelduseks oli, et ennustusjäägid peavad olema normaaljaotuslikud. Tuleb aga välja, et taolist lineaarset mudelit on võimalik üldistada ning eeldada seejuures väga erinevate jaotustega ennustusjääke (ennustusjäägid peavad olema eksponentide perekonnast). See võimaldab ühe tüüpmudeli abil modelleerida näiteks pidevtunnustega andmeid (normaaljaotus), loendusandmeid (Poissooni jaotus), eksponentseoseid, kui ka binaarseid andmeid (Bernoulli jaotus). Mudeli funktsionaalne vorm on järgnev:
\(E(y)=g^{-1}(X\beta)\)
kus:
Antud juhul keskendume jah/ei jaotusele ehk Bernoulli jaotusele, kus
\(g(\mu) = X\beta = ln\big(\frac{\mu}{1-\mu}\big)\) ja
\(\mu = \frac{1}{1+exp(-X\beta)}\)
Logistiline regressioon on olemuselt lineaarne klassifikaator, st ta sobitab kahe klassi vahele lineaarse eralduspiiri. Selle visualiseerimiseks võib teha simulatsiooni, kus võrdleme logistilist regressiooni mõne tuntuma mitte-lineaarse klassifikaatoriga (juhumets ehk random forest ja tugivektormasin ehk support vector machine).
library(randomForest)
library(e1071)
cos_func = function(x,y){
set.seed(200212)
x > y + cos(x)*rnorm(1,0,2)
}
pr = function(x) predict(x,test) %>% as.numeric() -1
data = expand.grid(x=seq(0,10,.25),
y=seq(0,10,.25)) %>%
mutate(label = cos_func(x,y))
log_fit = glm(label~.,data=data,family=binomial,control=glm.control(maxit=20))
rf_fit = randomForest(label~.,data=data,mtry=2)
svm_fit = svm(factor(label)~.,data=data,kernel="radial",probability=T,cost=5)
data %>%
mutate(logit_preds = predict(log_fit,data,"response")>.5,
rf_preds = predict(rf_fit,data,"response")>.5,
svm_preds = predict(svm_fit,data)) %>%
reshape2::melt(c("x","y")) %>%
ggplot() +
geom_point(aes(x,y,color=value)) +
coord_fixed() +
facet_wrap(~variable) +
theme_bw()
Valimisosalust mõõdetakse indiviiditasandil enamasti binaarsel kujul - kas vastaja hääletas või mitte. Kui soovime seda modelleerida, siis saame arvutada küsitlusele vastaja tõenäosuse (0-1) kuuluda kas ühte või teise rühma.
Valimisosaluse analüüsiks võtame ette 2015. aasta Riigikogu valimiste küsitluse ning püüame ennustada inimeste valimisosalust sotsiaal-demograafiliste tunnuste abil. Järjestame esmalt ordinaalsed tunnused ja tutvume andmestikuga.
palgad = c("Alla 200 €","200–275 €","276–350 €","351–425 €","426–500 €","501–575 €",
"576–700 €","701–1000 €","1001–1300 €","1301–1600 €","1601–1900 €",
"1901–2200 €","2201–2500 €","Üle 2500 €")
haridus = c("Põhiharidus","Keskharidus","Kõrgharidus")
RK_2015$wage = factor(RK_2015$wage,levels=palgad)
RK_2015$education = factor(RK_2015$education,levels=haridus)
summary(RK_2015)
vote age gender education
Ei hääletanud:150 Min. :18.00 Mees :350 Põhiharidus:118
Hääletas :646 1st Qu.:37.00 Naine:446 Keskharidus:417
Median :53.00 Kõrgharidus:261
Mean :51.97
3rd Qu.:67.00
Max. :90.00
ethnicity wage
Eestlane :653 701–1000 € :143
Mitte-eestlane:143 576–700 € :123
276–350 € : 96
1301–1600 €: 80
1001–1300 €: 78
351–425 € : 72
(Other) :204
Sõltuvtunnus on antud hetkel vote, mille puhul näeme, et 646 küsitletutest hääletas ja 150 ei hääletanud. Viime läbi logistilise regressiooni kõikide andmestikus kaasatud tunnustega. Logistilise regressiooni jaoks kasutame funktsiooni glm(...), kus tuleb lisaks valemile ja andmestikule märkida ennustusjääkide perekond family=binominal.
f = formula(vote ~ age + I(age^2) + gender + education + ethnicity + as.numeric(wage))
osalus_mudel = glm(f, data= RK_2015, family=binomial)
summary(osalus_mudel)
Call:
glm(formula = f, family = binomial, data = RK_2015)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.3317 0.3497 0.5297 0.6659 1.4727
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.537e+00 7.877e-01 -1.951 0.051055 .
age 3.143e-02 3.188e-02 0.986 0.324268
I(age^2) -2.829e-05 3.184e-04 -0.089 0.929201
genderNaine -8.873e-02 1.950e-01 -0.455 0.649071
educationKeskharidus 8.422e-01 2.441e-01 3.449 0.000562 ***
educationKõrgharidus 1.626e+00 3.057e-01 5.318 1.05e-07 ***
ethnicityMitte-eestlane 5.732e-03 2.471e-01 0.023 0.981495
as.numeric(wage) 9.712e-02 3.544e-02 2.740 0.006141 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 770.46 on 795 degrees of freedom
Residual deviance: 710.26 on 788 degrees of freedom
AIC: 726.26
Number of Fisher Scoring iterations: 5
Kuidas visualiseerida interaktsioonide marginaalefektie (ehk tingimuslik tõenäosus)? Sellist võimalust pakub funktsioon interplot(...) paketist interplot.
library(interplot)
interplot(osalus_mudel, var1 = "age", var2 = "age") +
theme_bw()
Rohkem võimalusi interploti kohta leiab siit: interplot: Plot the Effects of Variables in Interaction Terms
Et saada kätte lihtsamini tõlgendatavad šansside suhted, tuleb \(\beta\) väärtused eksponentida. \(\beta\) väärtused saame funktsiooniga coef() ning eksponendid exp(...).
(or = exp(coef(osalus_mudel)))
(Intercept) age I(age^2)
0.2150577 1.0319259 0.9999717
genderNaine educationKeskharidus educationKõrgharidus
0.9150965 2.3213569 5.0810969
ethnicityMitte-eestlane as.numeric(wage)
1.0057483 1.1019953
Šansside suhet ei saa siiski lineaarselt tõlgendada, sest konkreetne šansi suurus sõltub X väärtustest. Üks viis, kuidas kõikide X peale keskmist mõju suurust väljendada, on leida keskmine marginaalefekt. Selleks tuleb leida esmalt keskmine logaritm-suhte tihedus ning korrutada see koefitsientide suurustega. Kui see tundub praegu keeruline, siis antud hetkel ei olegi tarvis seda mõista. Huvi korral võib alustada tõenäosuse ja logaritm-suhte uurimist siit: Binary Outcome GLM Plots.
(ame = mean(dlogis(predict(osalus_mudel, type="link"))) * coef(osalus_mudel))
(Intercept) age I(age^2)
-2.163208e-01 4.423521e-03 -3.981973e-06
genderNaine educationKeskharidus educationKõrgharidus
-1.248868e-02 1.185380e-01 2.288028e-01
ethnicityMitte-eestlane as.numeric(wage)
8.067866e-04 1.367057e-02
Kui šansside suhet võis mõista “kui palju kasvab sündmuse šanss (\(\frac{y=1}{y=0}\)) X-i kasvades ühe võrra” (mitte segi ajada riskisuhtega (\(\frac{y=1}{y}\))!), siis keskmine marginaalefekt tähistab regressiooni joone keskmist kallakut.
Kuvame kõik leitud väärtused ühes tabelis kasutades selleks stargazerit.
stargazer(osalus_mudel,osalus_mudel,osalus_mudel,
dep.var.labels = "Hääletas valimistel",
column.labels = c("Beta koefitsiendid",
"Šansside suhe",
"Keskmine marginaalefekt"),
dep.var.caption = "",
coef = list(NULL, or, ame),
type="html")
| Hääletas valimistel | |||
| Beta koefitsiendid | Šansside suhe | Keskmine marginaalefekt | |
| (1) | (2) | (3) | |
| age | 0.031 | 1.032*** | 0.004 |
| (0.032) | (0.032) | (0.032) | |
| I(age2) | -0.00003 | 1.000*** | -0.00000 |
| (0.0003) | (0.0003) | (0.0003) | |
| genderNaine | -0.089 | 0.915*** | -0.012 |
| (0.195) | (0.195) | (0.195) | |
| educationKeskharidus | 0.842*** | 2.321*** | 0.119 |
| (0.244) | (0.244) | (0.244) | |
| educationKõrgharidus | 1.626*** | 5.081*** | 0.229 |
| (0.306) | (0.306) | (0.306) | |
| ethnicityMitte-eestlane | 0.006 | 1.006*** | 0.001 |
| (0.247) | (0.247) | (0.247) | |
| as.numeric(wage) | 0.097*** | 1.102*** | 0.014 |
| (0.035) | (0.035) | (0.035) | |
| Constant | -1.537* | 0.215 | -0.216 |
| (0.788) | (0.788) | (0.788) | |
| Observations | 796 | 796 | 796 |
| Log Likelihood | -355.130 | -355.130 | -355.130 |
| Akaike Inf. Crit. | 726.259 | 726.259 | 726.259 |
| Note: | p<0.1; p<0.05; p<0.01 | ||
Arvutame viimaks McFaddeni pseudo R-ruudu suuruse. Arvestada tuleks, et erinevalt lineaarse regressiooni R-ruudust, on antud pseudo R-ruudu suurus alati madalam (0.2-0.4 on juba hea tulemus).
\(R^2=1-\frac{logLik_{mod}}{logLik_{null}}\)
1 - osalus_mudel$deviance / osalus_mudel$null.deviance
[1] 0.07813415
Meenutame, et andmestikus on hääletanuid oluliselt rohkem kui mittehääletanuid. Mudelit saaks seega oluliselt parandada tasakaalustatud andmestikuga. Kasutame selleks paketti ROSE ja funktsiooni ovun.sample(...).
library(ROSE)
over = ovun.sample(vote~.,data=RK_2015,method="over",N=1292)$data
under = ovun.sample(vote~.,data=RK_2015,method="under",N=300)$data
both = ovun.sample(vote~.,data=RK_2015,method="both",p=.5,N=nrow(RK_2015))$data
osalus_over = glm(f, data=over, family=binomial)
osalus_under = glm(f, data=under, family=binomial)
osalus_both = glm(f, data=both, family=binomial)
1 - osalus_over$deviance / osalus_over$null.deviance
[1] 0.08366703
1 - osalus_under$deviance / osalus_under$null.deviance
[1] 0.07123121
1 - osalus_both$deviance / osalus_both$null.deviance
[1] 0.1038426
Hea. Pseudo R-ruut tõusis tasakaalustatud andmetega mõne protsendi võrra.
\(\chi^2\) test kontrollib, kas testidu mudel on oluliselt parem tühjast mudelist.
H&L test kontrollib, kas mudel sobitub halvasti. Hea on seega kõrge p-väärtus, mis näitab, et mudel sobitub hästi.
hoslem.test(osalus_mudel$y, fitted(osalus_mudel), g=10)
Hosmer and Lemeshow goodness of fit (GOF) test
data: osalus_mudel$y, fitted(osalus_mudel)
X-squared = 8.6897, df = 8, p-value = 0.3691
logit_preds = ifelse(predict(osalus_mudel,type="response")>.5,"Hääletas","Ei hääletanud")
confusionMatrix(logit_preds,RK_2015$vote)
Confusion Matrix and Statistics
Reference
Prediction Ei hääletanud Hääletas
Ei hääletanud 9 7
Hääletas 141 639
Accuracy : 0.8141
95% CI : (0.7853, 0.8405)
No Information Rate : 0.8116
P-Value [Acc > NIR] : 0.4496
Kappa : 0.0748
Mcnemar's Test P-Value : <2e-16
Sensitivity : 0.06000
Specificity : 0.98916
Pos Pred Value : 0.56250
Neg Pred Value : 0.81923
Prevalence : 0.18844
Detection Rate : 0.01131
Detection Prevalence : 0.02010
Balanced Accuracy : 0.52458
'Positive' Class : Ei hääletanud
ROC_curve = function(probs,response,cut=.5){
y_values = unique(response)
cutoffs = sort(unique(as.numeric(probs)))
tpr = sapply(cutoffs, function(cuts) sum((probs >= cuts) == T & response == y_values[1])/(sum(response == y_values[1])))
fpr = sapply(cutoffs, function(cuts) sum((probs >= cuts) == T & response == y_values[2])/(sum(response == y_values[2])))
stats = data.frame(cutoffs, tpr, fpr)
r = rep(y_values[2],length(response))
r[probs>cut] = y_values[1]
cm = confusionMatrix(r,response)
auc = cm$byClass["Balanced Accuracy"]
stats %>%
ggplot() +
geom_line(aes(fpr,tpr)) +
geom_ribbon(aes(x=fpr,ymax=tpr,ymin=fpr),alpha=.5,fill="steelblue") +
geom_label(aes(.5,.5,label=paste("AUC:",round(auc,3))))+
scale_x_continuous(breaks = seq(0,1,.1)) +
scale_y_continuous(breaks = seq(0,1,.1)) +
coord_fixed() +
theme_bw()
}
ROC_curve(predict(osalus_mudel,type="response"), RK_2015$vote)